q— title: “Experimental Data Analysis” author: “Savannah Weaver” date: “July 2021” output: pdf_document: toc: TRUE —
Packages
Data
This data was collected in the Spring of 2021 in conjuction with a study carried out in Cal Poly’s Herpetology class. Some lizards measured for that primary study were kept to observe physiological changes in response to different climate treatments.
Morphometrics & Hydration
Treatment Groups
variables: - individual lizard ID - temp_tmt_C = temperature treatment - humidity_tmt_percent = humidity treatment (high/low, not actually %) - trial_number = which set of lizards that individual was from - conclusion = how that individual’s experiment ended (died, canceled, or complete)
tmts <- read.csv("./data/exp_tmt_assignment.csv")Capture Data
variables: - date = date of capture & baseline measurements - individual lizard ID - mass_g = mass in grams - hematocrit_percent = % of blood sample that’s red blood cells - osmolality_mmol_kg = concentration of solutes in blood plasma - type = when the measurements were taken along the course of the experiment (all on capture day)
capture_hydration <- read.csv("./exported_data/capture_hydration.csv",
na.strings=c("","NA") # fix empty cells
) %>%
mutate(# correctly format date-only variable
date = as.Date(date, format = "%Y-%m-%d")
) %>%
# select only relevant variables
dplyr::select(date, individual_ID,
mass_g, hematocrit_percent, osmolality_mmol_kg
) %>%
dplyr::filter(individual_ID %in% tmts$individual_ID) %>%
mutate(type = as.factor("capture"))
summary(capture_hydration)## date individual_ID mass_g hematocrit_percent
## Min. :2021-04-19 Min. : 31.00 Min. : 8.20 Min. :16.00
## 1st Qu.:2021-04-26 1st Qu.: 57.00 1st Qu.:11.10 1st Qu.:32.75
## Median :2021-04-26 Median : 78.00 Median :12.65 Median :36.00
## Mean :2021-04-29 Mean : 77.46 Mean :12.18 Mean :36.08
## 3rd Qu.:2021-05-03 3rd Qu.: 98.25 3rd Qu.:13.32 3rd Qu.:39.00
## Max. :2021-05-10 Max. :122.00 Max. :15.00 Max. :54.00
## osmolality_mmol_kg type
## Min. :319.0 capture:52
## 1st Qu.:354.2
## Median :373.0
## Mean :373.7
## 3rd Qu.:392.2
## Max. :423.0extract SVL data separately from capture data:
SVL <- read.csv("./exported_data/capture_hydration.csv",
na.strings=c("","NA") # fix empty cells
) %>%
dplyr::select(individual_ID, SVL_mm) %>%
dplyr::filter(individual_ID %in% tmts$individual_ID)
summary(SVL)## individual_ID SVL_mm
## Min. : 31.00 Min. :59.00
## 1st Qu.: 57.00 1st Qu.:65.00
## Median : 78.00 Median :68.00
## Mean : 77.46 Mean :67.62
## 3rd Qu.: 98.25 3rd Qu.:70.00
## Max. :122.00 Max. :73.00extract capture CEWL cloacal temperature separately:
cap_CT <- read.csv("./exported_data/capture_hydration.csv",
na.strings=c("","NA") # fix empty cells
) %>%
dplyr::select(individual_ID, cloacal_temp_C) %>%
dplyr::filter(individual_ID %in% tmts$individual_ID)
summary(cap_CT)## individual_ID cloacal_temp_C
## Min. : 31.00 Min. :20.00
## 1st Qu.: 57.00 1st Qu.:22.00
## Median : 78.00 Median :24.00
## Mean : 77.46 Mean :23.68
## 3rd Qu.: 98.25 3rd Qu.:25.00
## Max. :122.00 Max. :28.00
## NA's :2Experiment Data
variables: - date = date of measurements - individual lizard ID - mass_g = mass in grams - hematocrit_percent = % of blood sample that’s red blood cells - osmolality_mmol_kg = concentration of solutes in blood plasma (mean of 1-3 replicates) - type = when the measurements were taken along the course of the experiment (either during experimental treatment or after rehab)
exp_dat <- read.csv("./data/experimental_data.csv",
na.strings=c("","NA") # fix empty cells
) %>%
# format date
dplyr::mutate(date = as.Date(date, format = "%m/%d/%y"),
type = as.factor(type)
) %>%
# select only variables to be analyzed
dplyr::select(date, individual_ID, mass_g,
hematocrit_percent, type,
osmolality_mmol_kg = osmolality_mmol_kg_replicate_mean)
summary(exp_dat)## date individual_ID mass_g hematocrit_percent
## Min. :2021-04-21 Min. : 31.00 Min. : 6.700 Min. :12.0
## 1st Qu.:2021-04-28 1st Qu.: 51.25 1st Qu.: 9.875 1st Qu.:23.0
## Median :2021-05-07 Median : 87.50 Median :11.250 Median :28.0
## Mean :2021-05-06 Mean : 77.85 Mean :11.076 Mean :27.8
## 3rd Qu.:2021-05-14 3rd Qu.:101.25 3rd Qu.:12.225 3rd Qu.:33.0
## Max. :2021-05-20 Max. :122.00 Max. :14.700 Max. :43.0
## NA's :19
## type osmolality_mmol_kg
## exp :98 Min. :298.0
## rehab:34 1st Qu.:342.0
## Median :355.0
## Mean :360.1
## 3rd Qu.:374.8
## Max. :441.0
## NA's :22Join Dataframes
Now, attach all the dataframes, only use individuals whose treatment was completed, and add a “day” variable for what day of treatment each lizard/observation was on. I also calculate SMI using the equation created in capture_analysis.
all_dat <- exp_dat %>%
# join data
rbind(capture_hydration) %>%
# add tmt group info
left_join(tmts, by = "individual_ID") %>%
dplyr::select(-notes) %>%
# add SVL value for each obs of each indiv.
# for computing BCI and scaled mass indices
left_join(SVL, by = "individual_ID") %>%
# only use completed experiment runs
dplyr::filter(conclusion == "complete") %>%
group_by(individual_ID) %>%
# reformat a lot of variables
mutate(capture_date = min(date),
day = as.numeric(date - capture_date),
humidity_tmt_percent = as.factor(humidity_tmt_percent),
individual_ID = as.factor(individual_ID),
temp_tmt_C = as.factor(temp_tmt_C),
trial_number = as.factor(trial_number),
conclusion = as.factor(conclusion),
SMI = mass_g * ((65.02158/SVL_mm) ^ (3.09059/sqrt(0.8944)))
)
summary(all_dat)## date individual_ID mass_g hematocrit_percent
## Min. :2021-04-19 37 : 6 Min. : 6.70 Min. :12.00
## 1st Qu.:2021-04-30 39 : 6 1st Qu.:10.20 1st Qu.:24.00
## Median :2021-05-07 40 : 6 Median :11.50 Median :30.00
## Mean :2021-05-06 49 : 6 Mean :11.27 Mean :29.58
## 3rd Qu.:2021-05-13 52 : 6 3rd Qu.:12.60 3rd Qu.:35.00
## Max. :2021-05-20 47 : 5 Max. :15.00 Max. :54.00
## (Other):116 NA's :12
## type osmolality_mmol_kg temp_tmt_C humidity_tmt_percent trial_number
## exp :82 Min. :298.0 25:151 dry :74 1:35
## rehab :34 1st Qu.:342.8 humid:77 2:24
## capture:35 Median :359.0 3:44
## Mean :362.5 4:48
## 3rd Qu.:379.0
## Max. :441.0
## NA's :15
## conclusion SVL_mm capture_date day
## complete:151 Min. :59.00 Min. :2021-04-19 Min. : 0.000
## 1st Qu.:66.00 1st Qu.:2021-04-26 1st Qu.: 2.000
## Median :68.00 Median :2021-05-03 Median : 4.000
## Mean :67.45 Mean :2021-04-30 Mean : 5.424
## 3rd Qu.:70.00 3rd Qu.:2021-05-10 3rd Qu.: 9.000
## Max. :73.00 Max. :2021-05-10 Max. :11.000
##
## SMI
## Min. : 7.343
## 1st Qu.: 8.990
## Median :10.011
## Mean : 9.983
## 3rd Qu.:10.751
## Max. :13.970
## unique(all_dat$individual_ID)## [1] 40 52 37 39 47 49 80 66 54 61 74 73 92 91 95 88 93 96 98
## [20] 89 99 81 97 104 108 122 118 109 113 105 114 101 117 102 103
## 35 Levels: 37 39 40 47 49 52 54 61 66 73 74 80 81 88 89 91 92 93 95 96 ... 122re-order some factors:
all_dat$humidity_tmt_percent <- factor(all_dat$humidity_tmt_percent,
levels = c("humid", "dry"),
labels = c("Humid", "Dry"))make a sub-dataframe without rehab data to prevent any mix-ups:
all_dat_no_rehab <- all_dat %>%
dplyr::filter(type != "rehab")Checks
Dates:
# check that capture dates are valid
unique(all_dat$capture_date)## [1] "2021-04-19" "2021-04-26" "2021-05-03" "2021-05-10"Check that each lizard only has an accurate number of measurements.
all_dat %>%
group_by(individual_ID, type) %>%
summarise(n = n()) %>%
arrange(type)## `summarise()` regrouping output by 'individual_ID' (override with `.groups` argument)## # A tibble: 104 x 3
## # Groups: individual_ID [35]
## individual_ID type n
## <fct> <fct> <int>
## 1 37 exp 4
## 2 39 exp 4
## 3 40 exp 4
## 4 47 exp 4
## 5 49 exp 4
## 6 52 exp 4
## 7 54 exp 2
## 8 61 exp 2
## 9 66 exp 2
## 10 73 exp 2
## # … with 94 more rowsThat all looks good, experimental measurements are either 4 (first trial) or 2 (other trials). I am excluding lizards that died in treatment from the analysis.
CEWL
Capture CEWL
variables: - date = date of capture & baseline measurements - individual lizard ID - region = which body area the measurement was taken from - TEWL_g_m2h = evaporative water loss - cloacal_temp_C = taken at measurement; influences CEWL
cap_CEWL <- read.csv("./exported_data/capture_CEWL.csv") %>%
dplyr::select(date, individual_ID, region, TEWL_g_m2h) %>%
mutate(#individual_ID = as.factor(individual_ID), # do later
date = as.Date(date, format = "%Y-%m-%d"),
region = as.factor(region),
day = as.factor("before"),
n_day = 0
) %>%
dplyr::filter(individual_ID %in% all_dat$individual_ID) %>%
left_join(cap_CT, by = 'individual_ID')
summary(cap_CEWL)## date individual_ID region TEWL_g_m2h day
## Min. :2021-04-19 Min. : 37.00 dewl:32 Min. : 7.48 before:163
## 1st Qu.:2021-04-26 1st Qu.: 73.00 dors:33 1st Qu.:20.54
## Median :2021-05-03 Median : 95.00 head:33 Median :27.43
## Mean :2021-05-02 Mean : 87.46 mite:32 Mean :29.30
## 3rd Qu.:2021-05-10 3rd Qu.:104.00 vent:33 3rd Qu.:36.91
## Max. :2021-05-10 Max. :122.00 Max. :62.94
## n_day cloacal_temp_C
## Min. :0 Min. :20.00
## 1st Qu.:0 1st Qu.:22.00
## Median :0 Median :24.00
## Mean :0 Mean :23.84
## 3rd Qu.:0 3rd Qu.:25.00
## Max. :0 Max. :28.00Post-Experiment CEWL
In the future, I could automate this like I did for the HOBO data.
Load in each of the post-rehab datafiles:
# trial 1
CEWL_t1 <- read.csv("./data/post_exp_CEWL/4-28-21-CEWL.csv", # filename
na.strings=c("","NA")) %>% # fix empty cells
# rename and select the pertinent variables/cols
# I have to do this for each one
# so they all have the same number of columns for joining
dplyr::select(date = Date,
Status,
ID = Comments,
TEWL_g_m2h = TEWL..g..m2h.. # rename
)
# trial 2
CEWL_t2 <- read.csv("./data/post_exp_CEWL/5-4-21-CEWL.csv",
na.strings=c("","NA")) %>%
dplyr::select(date = Date,
Status,
ID = Comments,
TEWL_g_m2h = TEWL..g..m2h..
)
# trial 3
CEWL_t3 <- read.csv("./data/post_exp_CEWL/5-11-21-CEWL.csv",
na.strings=c("","NA")) %>%
dplyr::select(date = Date,
Status,
ID = Comments,
TEWL_g_m2h = TEWL..g..m2h..
)
# trial 4
CEWL_t4 <- read.csv("./data/post_exp_CEWL/5-18-21-CEWL.csv",
na.strings=c("","NA")) %>%
dplyr::select(date = Date,
Status,
ID = Comments,
TEWL_g_m2h = TEWL..g..m2h..
)Load in cloacal temperatures:
exp_CT <- read.csv("./data/post_exp_CEWL_cloacal_temps.csv") %>%
mutate(date = as.Date(date, format = "%Y/%m/%d")) %>%
dplyr::select(-time)
summary(exp_CT)## date individual_ID cloacal_temp_C
## Min. :2021-04-28 Min. : 37.00 Min. :19.0
## 1st Qu.:2021-05-04 1st Qu.: 69.50 1st Qu.:21.0
## Median :2021-05-11 Median : 93.00 Median :23.0
## Mean :2021-05-09 Mean : 85.91 Mean :22.4
## 3rd Qu.:2021-05-18 3rd Qu.:103.50 3rd Qu.:23.0
## Max. :2021-05-18 Max. :122.00 Max. :26.0Join Dataframes
define not-in function:
`%nin%` = Negate(`%in%`)Merge all post-experiment CEWL, add cloacal temperature, add capture CEWL:
# merge all CEWL datafiles & reformat
CEWL <- CEWL_t1 %>% # trial 1
rbind(., CEWL_t2, # trial 2
CEWL_t3, # trial 3
CEWL_t4 # trial 4
) %>%
# remove any unsuccessful measurements
dplyr::filter(Status == "Normal") %>%
# extract individual_ID and region separately from the "ID" variable
separate(ID, c("individual_ID", "region")) %>%
# reformat data
dplyr::mutate(# reformat date
date = as.Date(date, format = "%m/%d/%y"),
# format individual ID
individual_ID = as.integer(individual_ID),
# set body region as a factor variable after getting only the consistent characters due to typos
region = as.factor(substring(region, 1, 4)),
# add when measurement taken
day = as.factor("after"),
n_day = 1 # technically day 8/9, just to help with figures
) %>%
# remove cols not relevant to stats
dplyr::select(-Status) %>%
# remove any rows with missing values
# none actually needed to be removed
dplyr::filter(complete.cases(.)) %>%
# add cloacal temperatures
left_join(exp_CT, by = c("date", "individual_ID")) %>%
# now matching dataframes, add capture CEWL data
rbind(cap_CEWL) %>%
# add tmt assignments
left_join(tmts, by = "individual_ID") %>%
mutate(humidity_tmt_percent = as.factor(humidity_tmt_percent),
individual_ID = as.factor(individual_ID),
conclusion = as.factor(conclusion),
trial_number = as.factor(trial_number)
) %>%
# lizards 49 & 80 are missing pre-exp CEWL, so remove them
dplyr::filter((individual_ID %nin% c('49', '80')))
# every lizard should have 10 measurements
summary(CEWL)## date individual_ID region TEWL_g_m2h day
## Min. :2021-04-19 37 : 10 dewl:65 Min. : 4.60 after :163
## 1st Qu.:2021-05-03 39 : 10 dors:65 1st Qu.: 20.09 before:163
## Median :2021-05-10 40 : 10 head:66 Median : 27.18
## Mean :2021-05-06 47 : 10 mite:64 Mean : 30.69
## 3rd Qu.:2021-05-11 52 : 10 vent:66 3rd Qu.: 38.72
## Max. :2021-05-18 54 : 10 Max. :106.38
## (Other):266
## n_day cloacal_temp_C temp_tmt_C humidity_tmt_percent trial_number
## Min. :0.0 Min. :19.00 Min. :25 dry :158 1: 50
## 1st Qu.:0.0 1st Qu.:21.00 1st Qu.:25 humid:168 2: 48
## Median :0.5 Median :23.00 Median :25 3:110
## Mean :0.5 Mean :23.11 Mean :25 4:118
## 3rd Qu.:1.0 3rd Qu.:24.75 3rd Qu.:25
## Max. :1.0 Max. :28.00 Max. :25
##
## conclusion notes
## complete:326 Length:326
## Class :character
## Mode :character
##
##
##
## Check that data looks correct:
CEWL %>%
group_by(individual_ID, day) %>%
summarise(n = n()) %>%
arrange(individual_ID, n)## `summarise()` regrouping output by 'individual_ID' (override with `.groups` argument)## # A tibble: 66 x 3
## # Groups: individual_ID [33]
## individual_ID day n
## <fct> <fct> <int>
## 1 37 after 5
## 2 37 before 5
## 3 39 after 5
## 4 39 before 5
## 5 40 after 5
## 6 40 before 5
## 7 47 after 5
## 8 47 before 5
## 9 52 after 5
## 10 52 before 5
## # … with 56 more rowsEverything looks great! (after removing the observations for the two lizards with missing pre-experiment CEWL measurements.)
Before/after aren’t perfectly even because sometimes we were unable to get the AquaFlux to equilibrate and take a measurement.
Finally, make a small edit so the regions are spelled out completely. This requires reordering factor levels:
CEWL$region <- factor(CEWL$region,
levels = c("dors", "vent", "head", "dewl", "mite"),
labels = c("Dorsum", "Ventrum", "Head",
"Dewlap", "Mite Patch")
)
CEWL$humidity_tmt_percent <- factor(CEWL$humidity_tmt_percent,
levels = c("humid", "dry"),
labels = c("Humid", "Dry"))
CEWL$day <- factor(CEWL$day,
levels = c("before", "after"),
labels = c("Before", "After"))
summary(CEWL)## date individual_ID region TEWL_g_m2h
## Min. :2021-04-19 37 : 10 Dorsum :65 Min. : 4.60
## 1st Qu.:2021-05-03 39 : 10 Ventrum :66 1st Qu.: 20.09
## Median :2021-05-10 40 : 10 Head :66 Median : 27.18
## Mean :2021-05-06 47 : 10 Dewlap :65 Mean : 30.69
## 3rd Qu.:2021-05-11 52 : 10 Mite Patch:64 3rd Qu.: 38.72
## Max. :2021-05-18 54 : 10 Max. :106.38
## (Other):266
## day n_day cloacal_temp_C temp_tmt_C humidity_tmt_percent
## Before:163 Min. :0.0 Min. :19.00 Min. :25 Humid:168
## After :163 1st Qu.:0.0 1st Qu.:21.00 1st Qu.:25 Dry :158
## Median :0.5 Median :23.00 Median :25
## Mean :0.5 Mean :23.11 Mean :25
## 3rd Qu.:1.0 3rd Qu.:24.75 3rd Qu.:25
## Max. :1.0 Max. :28.00 Max. :25
##
## trial_number conclusion notes
## 1: 50 complete:326 Length:326
## 2: 48 Class :character
## 3:110 Mode :character
## 4:118
##
##
## Export Data Frames for Power Analyses
#write.csv(all_dat, "exported_data/exp_effects_hydration.csv")
#write.csv(CEWL, "exported_data/exp_effects_CEWL.csv")Data Distributions
Histograms
Mass
all_dat %>%
ggplot(., aes(x = mass_g)) +
geom_histogram(color = "black", fill="steelblue", bins=50) +
theme_classic() +
xlab("Mass (g)") +
ylab("Count")
simple.eda(all_dat$mass_g)
shapiro.test(all_dat$mass_g)##
## Shapiro-Wilk normality test
##
## data: all_dat$mass_g
## W = 0.97747, p-value = 0.01396Mass distribution not normal, skewed to the left.
Scaled Mass Index
all_dat %>%
ggplot(., aes(x = SMI)) +
geom_histogram(color = "black", fill="steelblue", bins=50) +
theme_classic() +
xlab("Scaled Mass Index (g)") +
ylab("Count")
simple.eda(all_dat$SMI)
shapiro.test(all_dat$SMI)##
## Shapiro-Wilk normality test
##
## data: all_dat$SMI
## W = 0.99012, p-value = 0.3712Hematocrit
all_dat %>%
ggplot(., aes(x = hematocrit_percent)) +
geom_histogram(color = "black", fill="steelblue", bins=50) +
theme_classic() +
xlab("Hematocrit (%)") +
ylab("Count")## Warning: Removed 12 rows containing non-finite values (stat_bin).
simple.eda(all_dat$hematocrit_percent)
shapiro.test(all_dat$hematocrit_percent)##
## Shapiro-Wilk normality test
##
## data: all_dat$hematocrit_percent
## W = 0.98984, p-value = 0.4089Visually, looks slightly skewed to the right, but statistically, the distribution of hematocrit is normal.
Osmolality
all_dat %>%
ggplot(., aes(x = osmolality_mmol_kg)) +
geom_histogram(color = "black", fill="steelblue", bins=50) +
theme_classic() +
xlab("Osmolality (mmol / kg)") +
ylab("Count")## Warning: Removed 15 rows containing non-finite values (stat_bin).
simple.eda(all_dat$osmolality_mmol_kg)
shapiro.test(all_dat$osmolality_mmol_kg)##
## Shapiro-Wilk normality test
##
## data: all_dat$osmolality_mmol_kg
## W = 0.98331, p-value = 0.09544Visually, looks slightly skewed to the right, but statistically, the distribution of osmolality is normal.
Data Removal
test the effect of removing the day 2/6 data, for consistency’s sake
trimmed <- all_dat_no_rehab %>%
dplyr::filter(day %in% c(0, 4, 8, 9))
summary(trimmed)## date individual_ID mass_g hematocrit_percent
## Min. :2021-04-19 37 : 3 Min. : 6.70 Min. :12.00
## 1st Qu.:2021-04-30 39 : 3 1st Qu.:10.20 1st Qu.:27.00
## Median :2021-05-07 40 : 3 Median :11.60 Median :32.00
## Mean :2021-05-05 47 : 3 Mean :11.44 Mean :31.39
## 3rd Qu.:2021-05-11 49 : 3 3rd Qu.:12.80 3rd Qu.:36.00
## Max. :2021-05-18 52 : 3 Max. :15.00 Max. :54.00
## (Other):87
## type osmolality_mmol_kg temp_tmt_C humidity_tmt_percent trial_number
## exp :70 Min. :313.0 25:105 Humid:54 1:18
## rehab : 0 1st Qu.:341.0 Dry :51 2:18
## capture:35 Median :356.5 3:33
## Mean :361.0 4:36
## 3rd Qu.:376.8
## Max. :426.0
## NA's :3
## conclusion SVL_mm capture_date day
## complete:105 Min. :59.0 Min. :2021-04-19 Min. :0.000
## 1st Qu.:66.0 1st Qu.:2021-04-26 1st Qu.:0.000
## Median :68.0 Median :2021-05-03 Median :4.000
## Mean :67.6 Mean :2021-05-01 Mean :4.057
## 3rd Qu.:70.0 3rd Qu.:2021-05-10 3rd Qu.:8.000
## Max. :73.0 Max. :2021-05-10 Max. :9.000
##
## SMI
## Min. : 7.604
## 1st Qu.: 9.142
## Median :10.052
## Mean :10.061
## 3rd Qu.:10.884
## Max. :13.970
## re-order: NOTE: same difference to make categorical or continuous once day 2/6 are removed
trimmed$day <- factor(trimmed$day,
levels = c("0", "4", "8", "9"),
labels = c("Before", "Mid", "After", "After"))Basic Figures & Models
Figure Overlays
all_dat_trimmed_means <- trimmed %>%
group_by(humidity_tmt_percent, day) %>%
summarise(SMI_mean = mean(SMI),
hct_mean = mean(hematocrit_percent))## `summarise()` regrouping output by 'humidity_tmt_percent' (override with `.groups` argument)all_dat_trimmed_mean_osml <- trimmed %>%
dplyr::filter(complete.cases(osmolality_mmol_kg)) %>%
group_by(humidity_tmt_percent, day) %>%
summarise(osml_mean = mean(osmolality_mmol_kg))## `summarise()` regrouping output by 'humidity_tmt_percent' (override with `.groups` argument)Mass ~ Time
I won’t be using this, SMI is more applicable.
Just look at plot:
all_dat_no_rehab %>%
ggplot(data = .) +
geom_point(aes(x = day,
y = mass_g,
color = humidity_tmt_percent
),
size = 1,
alpha = 0.6) +
stat_smooth(aes(x = day,
y = mass_g,
color = humidity_tmt_percent
),
formula = y ~ x,
method = "lm",
se = F,
size = 1.6,
alpha = 1 ) +
theme_classic() +
scale_x_continuous(breaks = c(0, 2, 4, 6, 8)) +
scale_color_brewer(palette = "Set2",
name = "Humidity Treatment") +
xlab("Day") +
ylab("Mass (g)") +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 10),
legend.text.align = 0
)
SMI ~ Time
plot over course of experiment:
ggplot() +
geom_point(data = trimmed,
aes(x = day,
y = SMI,
color = humidity_tmt_percent
),
size = 1,
alpha = 0.6) +
geom_line(data = trimmed,
aes(x = day,
y = SMI,
group = individual_ID,
color = humidity_tmt_percent),
alpha = 0.2) +
geom_line(data = all_dat_trimmed_means,
aes(x = day,
y = SMI_mean,
group = humidity_tmt_percent,
color = humidity_tmt_percent),
size = 1.6,
alpha = 1) +
theme_classic() +
scale_color_brewer(palette = "Set2",
name = "Humidity Treatment") +
xlab("") +
ylab("Body Condition (g)") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text = element_text(color = "black",
family = "sans",
size = 18),
legend.text.align = 0,
legend.position = "none"
) -> tmt_effects_SMI
tmt_effects_SMI
# export figure
ggsave(filename = "tmt_effects_SMI.jpeg",
plot = tmt_effects_SMI,
path = "./final_figures",
device = "jpeg",
dpi = 1200,
width = 5, height = 4)Check whether means started out different:
SMI_diff_lm <- all_dat_no_rehab %>%
dplyr::filter(day == 0) %>%
lm(data = ., SMI ~ humidity_tmt_percent)
summary(SMI_diff_lm)##
## Call:
## lm(formula = SMI ~ humidity_tmt_percent, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9405 -0.7429 -0.0401 0.7385 3.2183
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.7517 0.2811 38.25 <2e-16 ***
## humidity_tmt_percentDry -0.2904 0.4033 -0.72 0.476
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.192 on 33 degrees of freedom
## Multiple R-squared: 0.01547, Adjusted R-squared: -0.01436
## F-statistic: 0.5187 on 1 and 33 DF, p-value: 0.4765NOT significantly different, which is good.
Check whether means ended differently:
SMI_diff_lm_end <- all_dat_no_rehab %>%
dplyr::filter(day %in% c(8, 9)) %>%
lm(data = ., SMI ~ humidity_tmt_percent)
summary(SMI_diff_lm_end)##
## Call:
## lm(formula = SMI ~ humidity_tmt_percent, data = .)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.85144 -0.69873 -0.07453 0.80895 2.31088
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.8852 0.2585 38.246 <2e-16 ***
## humidity_tmt_percentDry -0.6287 0.3709 -1.695 0.0994 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.097 on 33 degrees of freedom
## Multiple R-squared: 0.08012, Adjusted R-squared: 0.05224
## F-statistic: 2.874 on 1 and 33 DF, p-value: 0.09943mixed model:
SMI_mod <- lme4::lmer(data = all_dat_no_rehab,
SMI ~ day*humidity_tmt_percent +
(1|trial_number))
summary(SMI_mod)## Linear mixed model fit by REML ['lmerMod']
## Formula: SMI ~ day * humidity_tmt_percent + (1 | trial_number)
## Data: all_dat_no_rehab
##
## REML criterion at convergence: 369.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.84407 -0.68219 -0.04622 0.62984 2.74822
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 0.01241 0.1114
## Residual 1.27504 1.1292
## Number of obs: 117, groups: trial_number, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 10.79499 0.24058 44.871
## day -0.10598 0.04510 -2.350
## humidity_tmt_percentDry -0.40660 0.33497 -1.214
## day:humidity_tmt_percentDry -0.04498 0.06463 -0.696
##
## Correlation of Fixed Effects:
## (Intr) day hmd__D
## day -0.759
## hmdty_tmt_D -0.678 0.545
## dy:hmdty__D 0.530 -0.698 -0.782drop1(SMI_mod) # this was more intuitive using lme4::lmer but I switched to be able to get p-values## Single term deletions
##
## Model:
## SMI ~ day * humidity_tmt_percent + (1 | trial_number)
## npar AIC
## <none> 369.24
## day:humidity_tmt_percent 1 367.74# drop interaction term
SMI_mod2 <- lmerTest::lmer(data = all_dat_no_rehab,
SMI ~ day + humidity_tmt_percent +
(1|trial_number))
summary(SMI_mod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SMI ~ day + humidity_tmt_percent + (1 | trial_number)
## Data: all_dat_no_rehab
##
## REML criterion at convergence: 366.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.84939 -0.64254 -0.07319 0.69627 2.67483
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 0.01262 0.1123
## Residual 1.26911 1.1265
## Number of obs: 117, groups: trial_number, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.88366 0.20374 20.15193 53.420 < 2e-16 ***
## day -0.12788 0.03223 110.92424 -3.967 0.000129 ***
## humidity_tmt_percentDry -0.58888 0.20841 111.10846 -2.826 0.005597 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) day
## day -0.641
## hmdty_tmt_D -0.498 0.000#drop1(SMI_mod2)
SMI_mod2t <- lmerTest::lmer(data = trimmed,
SMI ~ day + humidity_tmt_percent +
(1|trial_number))## boundary (singular) fit: see ?isSingularsummary(SMI_mod2t)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: SMI ~ day + humidity_tmt_percent + (1 | trial_number)
## Data: trimmed
##
## REML criterion at convergence: 327.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.76998 -0.68801 -0.00593 0.73878 2.75012
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 0.000 0.000
## Residual 1.303 1.141
## Number of obs: 105, groups: trial_number, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 10.8313 0.2212 101.0000 48.964 < 2e-16 ***
## dayMid -0.6173 0.2728 101.0000 -2.263 0.025803 *
## dayAfter -1.0308 0.2728 101.0000 -3.778 0.000267 ***
## humidity_tmt_percentDry -0.4543 0.2229 101.0000 -2.039 0.044091 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dayMid dyAftr
## dayMid -0.617
## dayAfter -0.617 0.500
## hmdty_tmt_D -0.489 0.000 0.000
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular#drop1(SMI_mod2)SMI is best predicted by day and treatment, but not including their interaction.
write.csv(broom.mixed::tidy(SMI_mod2t),
"./best models/exp_effects_SMI.csv")Hct ~ Time
ggplot() +
geom_point(data = trimmed,
aes(x = day,
y = hematocrit_percent,
color = humidity_tmt_percent
),
size = 1,
alpha = 0.6) +
geom_line(data = trimmed,
aes(x = day,
y = hematocrit_percent,
group = individual_ID,
color = humidity_tmt_percent),
alpha = 0.2) +
geom_line(data = all_dat_trimmed_means,
aes(x = day,
y = hct_mean,
group = humidity_tmt_percent,
color = humidity_tmt_percent),
size = 1.6,
alpha = 1) +
theme_classic() +
scale_color_brewer(palette = "Set2",
name = "Humidity Treatment") +
xlab("Day") +
ylab("Hematocrit (%)") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text.align = 0,
legend.position = "none"
) -> tmt_effects_hct
tmt_effects_hct
# export figure
ggsave(filename = "tmt_effects_hct.jpeg",
plot = tmt_effects_hct,
path = "./final_figures",
device = "jpeg",
dpi = 1200,
width = 5, height = 4)this model seemed to work well with indiv as a random factor, but still excluded because it’s probably unnecessary
hct_mod <- all_dat_no_rehab %>%
dplyr::filter(complete.cases(hematocrit_percent)) %>%
lme4::lmer(data = .,
hematocrit_percent ~ day + humidity_tmt_percent +
(1|trial_number))
summary(hct_mod)## Linear mixed model fit by REML ['lmerMod']
## Formula: hematocrit_percent ~ day + humidity_tmt_percent + (1 | trial_number)
## Data: .
##
## REML criterion at convergence: 765.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2154 -0.6104 0.0919 0.6453 2.7070
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 8.339 2.888
## Residual 40.121 6.334
## Number of obs: 117, groups: trial_number, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 35.2145 1.8208 19.340
## day -1.0534 0.1813 -5.811
## humidity_tmt_percentDry -0.2733 1.1724 -0.233
##
## Correlation of Fixed Effects:
## (Intr) day
## day -0.403
## hmdty_tmt_D -0.314 0.000drop1(hct_mod)## Single term deletions
##
## Model:
## hematocrit_percent ~ day + humidity_tmt_percent + (1 | trial_number)
## npar AIC
## <none> 778.55
## day 1 806.68
## humidity_tmt_percent 1 776.61# drop humidity
hct_mod2 <- all_dat_no_rehab %>%
dplyr::filter(complete.cases(hematocrit_percent)) %>%
lmerTest::lmer(data = .,
hematocrit_percent ~ day +
(1|trial_number))
summary(hct_mod2)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: hematocrit_percent ~ day + (1 | trial_number)
## Data: .
##
## REML criterion at convergence: 767.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.2511 -0.5977 0.1135 0.6655 2.7398
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 8.36 2.891
## Residual 39.78 6.307
## Number of obs: 117, groups: trial_number, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 35.0810 1.7278 4.5115 20.304 1.31e-05 ***
## day -1.0534 0.1805 112.0809 -5.836 5.29e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## day -0.423drop1(hct_mod2)## Single term deletions using Satterthwaite's method:
##
## Model:
## hematocrit_percent ~ day + (1 | trial_number)
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## day 1354.8 1354.8 1 112.08 34.057 5.287e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1hct_mod2t <- trimmed %>%
dplyr::filter(complete.cases(hematocrit_percent)) %>%
lmerTest::lmer(data = .,
hematocrit_percent ~ day +
(1|trial_number))
summary(hct_mod2t)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: hematocrit_percent ~ day + (1 | trial_number)
## Data: .
##
## REML criterion at convergence: 676.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.14009 -0.58850 0.00871 0.64072 2.74219
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 7.253 2.693
## Residual 37.932 6.159
## Number of obs: 105, groups: trial_number, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 35.592 1.712 4.853 20.794 6.24e-06 ***
## dayMid -5.771 1.472 98.766 -3.920 0.000163 ***
## dayAfter -8.057 1.472 98.766 -5.473 3.36e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dayMid
## dayMid -0.430
## dayAfter -0.430 0.500The model AIC is slightly better without the interaction effect, so I removed that. The effect of humidity could ALSO be dropped, so humidity treatment was not an important factor affecting hematocrit, but how many days lizards were in treatment was. Both treatment groups lost hematocrit at approximately the same rate.
write.csv(broom::tidy(hct_mod2t),
"./best models/exp_effects_hct.csv")Osml ~ Time
ggplot() +
geom_point(data = trimmed,
aes(x = day,
y = osmolality_mmol_kg,
color = humidity_tmt_percent
),
size = 1,
alpha = 0.6) +
geom_line(data = trimmed,
aes(x = day,
y = osmolality_mmol_kg,
group = individual_ID,
color = humidity_tmt_percent),
alpha = 0.2) +
geom_line(data = all_dat_trimmed_mean_osml,
aes(x = day,
y = osml_mean,
group = humidity_tmt_percent,
color = humidity_tmt_percent),
size = 1.6,
alpha = 1) +
theme_classic() +
scale_color_brewer(palette = "Set2",
name = "Humidity Treatment") +
xlab("") +
ylab("Osmolality (mmol / kg)") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text.align = 0,
legend.position = "none"
) -> tmt_effects_osml
tmt_effects_osml## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).
# export figure
ggsave(filename = "tmt_effects_osml.jpeg",
plot = tmt_effects_osml,
path = "./final_figures",
device = "jpeg",
dpi = 1200,
width = 5, height = 4)## Warning: Removed 3 rows containing missing values (geom_point).
## Warning: Removed 3 row(s) containing missing values (geom_path).singular warning - do NOT include individual ID as a random effect
osml_mod <- all_dat_no_rehab %>%
dplyr::filter(complete.cases(osmolality_mmol_kg)) %>%
lmerTest::lmer(data = .,
osmolality_mmol_kg ~ day * humidity_tmt_percent +
(1|trial_number))
summary(osml_mod)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: osmolality_mmol_kg ~ day * humidity_tmt_percent + (1 | trial_number)
## Data: .
##
## REML criterion at convergence: 1018.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1763 -0.6722 -0.1769 0.6925 2.5285
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 302.6 17.40
## Residual 469.3 21.66
## Number of obs: 114, groups: trial_number, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 369.938 9.804 4.188 37.733 1.83e-06 ***
## day -1.315 0.877 106.917 -1.499 0.137
## humidity_tmt_percentDry -8.152 6.439 106.912 -1.266 0.208
## day:humidity_tmt_percentDry 1.901 1.266 106.908 1.502 0.136
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) day hmd__D
## day -0.357
## hmdty_tmt_D -0.319 0.542
## dy:hmdty__D 0.246 -0.691 -0.776drop1(osml_mod)## Single term deletions using Satterthwaite's method:
##
## Model:
## osmolality_mmol_kg ~ day * humidity_tmt_percent + (1 | trial_number)
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## day:humidity_tmt_percent 1058.1 1058.1 1 106.91 2.2545 0.1362osml_modt <- trimmed %>%
dplyr::filter(complete.cases(osmolality_mmol_kg)) %>%
lmerTest::lmer(data = .,
osmolality_mmol_kg ~ day * humidity_tmt_percent +
(1|trial_number))
summary(osml_modt)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: osmolality_mmol_kg ~ day * humidity_tmt_percent + (1 | trial_number)
## Data: .
##
## REML criterion at convergence: 872.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.66849 -0.79364 0.00361 0.51469 2.41228
##
## Random effects:
## Groups Name Variance Std.Dev.
## trial_number (Intercept) 313.5 17.71
## Residual 393.8 19.84
## Number of obs: 102, groups: trial_number, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 378.509 10.035 4.379 37.719 1.13e-06
## dayMid -29.500 6.615 92.951 -4.460 2.29e-05
## dayAfter -13.885 6.713 92.956 -2.068 0.0414
## humidity_tmt_percentDry -10.215 6.713 92.956 -1.522 0.1315
## dayMid:humidity_tmt_percentDry 12.324 9.491 92.951 1.298 0.1974
## dayAfter:humidity_tmt_percentDry 17.783 9.721 92.954 1.829 0.0705
##
## (Intercept) ***
## dayMid ***
## dayAfter *
## humidity_tmt_percentDry
## dayMid:humidity_tmt_percentDry
## dayAfter:humidity_tmt_percentDry .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dayMid dyAftr hmd__D dM:__D
## dayMid -0.330
## dayAfter -0.326 0.493
## hmdty_tmt_D -0.326 0.493 0.485
## dyMd:hmd__D 0.230 -0.697 -0.343 -0.707
## dyAftr:h__D 0.224 -0.340 -0.690 -0.690 0.488drop1(osml_modt)## Single term deletions using Satterthwaite's method:
##
## Model:
## osmolality_mmol_kg ~ day * humidity_tmt_percent + (1 | trial_number)
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## day:humidity_tmt_percent 1402.8 701.41 2 92.953 1.7812 0.1741The model seems good as-is.
write.csv(broom::tidy(osml_modt),
"./best models/exp_effects_osml.csv")Change in Osmolality
osml_d0 <- all_dat_no_rehab %>%
dplyr::filter(day == 0) %>%
dplyr::select(individual_ID, osml0 = osmolality_mmol_kg,
humidity_tmt_percent)
osml_d8 <- all_dat_no_rehab %>%
dplyr::filter(day %in% c(8,9)) %>%
dplyr::select(individual_ID, osml89 = osmolality_mmol_kg,
humidity_tmt_percent)
osml_diffs <- osml_d0 %>%
left_join(osml_d8) %>%
mutate(osml_change = osml89 - osml0)## Joining, by = c("individual_ID", "humidity_tmt_percent")boxplot:
osml_diffs %>%
ggplot(data = .) +
geom_boxplot(aes(x = humidity_tmt_percent,
y = osml_change,
group = humidity_tmt_percent,
color = humidity_tmt_percent
),
size = 1,
alpha = 1) +
theme_classic() +
geom_hline(yintercept = 0) +
xlab("") +
ylab("Osmolality Change (mmol / kg)") +
scale_color_brewer(palette = "Set2") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 14),
legend.text.align = 0,
legend.position = "none"
)## Warning: Removed 3 rows containing non-finite values (stat_boxplot).
CEWL ~ Before/After
try a boxplot:
CEWL %>%
ggplot(data = .) +
geom_boxplot(aes(x = humidity_tmt_percent,
y = TEWL_g_m2h,
color = day
),
size = 1,
alpha = 1) +
facet_wrap(~region) +
theme_classic() +
xlab("") +
ylab("CEWL (g / m^2h)") +
scale_color_brewer(palette = "Set2") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 14),
legend.text.align = 0,
legend.position = "bottom"
)
this is difficult to see changes, I think a line graph would be better…
CEWL %>%
ggplot(data = .) +
geom_point(aes(x = n_day,
y = TEWL_g_m2h,
color = humidity_tmt_percent
),
size = 1,
alpha = 0.6) +
geom_line(aes(x = n_day,
y = TEWL_g_m2h,
group = individual_ID,
color = humidity_tmt_percent),
alpha = 0.2) +
stat_smooth(aes(x = n_day,
y = TEWL_g_m2h,
color = humidity_tmt_percent
),
formula = y ~ x,
method = "lm",
se = F,
size = 1.6,
alpha = 1) +
theme_classic() +
scale_color_brewer(palette = "Set2",
name = "Humidity Treatment") +
facet_wrap(~region, ncol = 2) +
scale_x_continuous(breaks = c(0, 1),
labels = c("0" = "before", "1" = "after")
) +
xlab("") +
ylab(bquote('CEWL (g/'*m^2*'h)')) +
theme(text = element_text(color = "black",
family = "sans",
size = 12),
axis.text = element_text(color = "black",
family = "sans",
size = 10,
angle = 90),
legend.text = element_text(color = "black",
family = "sans",
size = 10),
legend.text.align = 0,
legend.position = c(0.75,0.12),
#legend.justification = c(1, 1)
) -> CEWL_tmt_fig
CEWL_tmt_fig
# export figure
#ggsave(filename = "tmt_effects_CEWL.jpeg",
# plot = CEWL_tmt_fig,
# path = "./final_figures",
# device = "jpeg",
# dpi = 1200,
# width = 4, height = 6)I saved the legend separately to make the figure layout better.
CEWL_mod <- CEWL %>%
lme4::lmer(data = .,
TEWL_g_m2h ~ day * humidity_tmt_percent * region +
cloacal_temp_C +
(1|trial_number/individual_ID))
summary(CEWL_mod)## Linear mixed model fit by REML ['lmerMod']
## Formula: TEWL_g_m2h ~ day * humidity_tmt_percent * region + cloacal_temp_C +
## (1 | trial_number/individual_ID)
## Data: .
##
## REML criterion at convergence: 2441.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4361 -0.5707 -0.0822 0.4555 4.1852
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID:trial_number (Intercept) 38.212 6.182
## trial_number (Intercept) 8.141 2.853
## Residual 123.274 11.103
## Number of obs: 326, groups: individual_ID:trial_number, 33; trial_number, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -69.91895 11.98454 -5.834
## dayAfter 20.55665 3.96057 5.190
## humidity_tmt_percentDry 0.09390 4.44521 0.021
## regionVentrum 14.06471 3.80826 3.693
## regionHead 5.20176 3.80826 1.366
## regionDewlap 2.76945 3.87231 0.715
## regionMite Patch 5.02118 3.80826 1.318
## cloacal_temp_C 3.90353 0.48902 7.982
## dayAfter:humidity_tmt_percentDry -22.17103 5.52280 -4.014
## dayAfter:regionVentrum 4.82016 5.43126 0.887
## dayAfter:regionHead -7.17102 5.43126 -1.320
## dayAfter:regionDewlap -2.93470 5.46212 -0.537
## dayAfter:regionMite Patch -6.04716 5.47797 -1.104
## humidity_tmt_percentDry:regionVentrum -0.34596 5.46920 -0.063
## humidity_tmt_percentDry:regionHead -0.03395 5.46920 -0.006
## humidity_tmt_percentDry:regionDewlap 0.90930 5.51399 0.165
## humidity_tmt_percentDry:regionMite Patch 2.90477 5.52066 0.526
## dayAfter:humidity_tmt_percentDry:regionVentrum -2.37516 7.76641 -0.306
## dayAfter:humidity_tmt_percentDry:regionHead 5.28070 7.76641 0.680
## dayAfter:humidity_tmt_percentDry:regionDewlap 5.66684 7.82338 0.724
## dayAfter:humidity_tmt_percentDry:regionMite Patch 1.47809 7.83657 0.189##
## Correlation matrix not shown by default, as p = 21 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itdrop1(CEWL_mod)## Single term deletions
##
## Model:
## TEWL_g_m2h ~ day * humidity_tmt_percent * region + cloacal_temp_C +
## (1 | trial_number/individual_ID)
## npar AIC
## <none> 2570.1
## cloacal_temp_C 1 2628.2
## day:humidity_tmt_percent:region 4 2563.8Drop triple interaction. I think the day:region standalone would be weird too.
CEWL_mod2 <- CEWL %>%
lme4::lmer(data = .,
TEWL_g_m2h ~
humidity_tmt_percent * (day + region) +
cloacal_temp_C +
(1|trial_number/individual_ID))
summary(CEWL_mod2)## Linear mixed model fit by REML ['lmerMod']
## Formula: TEWL_g_m2h ~ humidity_tmt_percent * (day + region) + cloacal_temp_C +
## (1 | trial_number/individual_ID)
## Data: .
##
## REML criterion at convergence: 2490.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.4790 -0.6180 -0.0854 0.4480 4.0055
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID:trial_number (Intercept) 38.017 6.166
## trial_number (Intercept) 8.677 2.946
## Residual 123.628 11.119
## Number of obs: 326, groups: individual_ID:trial_number, 33; trial_number, 4
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) -68.6125 11.9022 -5.765
## humidity_tmt_percentDry -0.8482 3.7242 -0.228
## dayAfter 18.2857 1.9294 9.478
## regionVentrum 16.5143 2.7189 6.074
## regionHead 1.6557 2.7189 0.609
## regionDewlap 1.3212 2.7199 0.486
## regionMite Patch 2.0988 2.7414 0.766
## cloacal_temp_C 3.8933 0.4894 7.956
## humidity_tmt_percentDry:dayAfter -20.1573 2.4933 -8.084
## humidity_tmt_percentDry:regionVentrum -1.5730 3.8883 -0.405
## humidity_tmt_percentDry:regionHead 2.5669 3.8883 0.660
## humidity_tmt_percentDry:regionDewlap 3.6718 3.9063 0.940
## humidity_tmt_percentDry:regionMite Patch 3.4656 3.9211 0.884##
## Correlation matrix not shown by default, as p = 13 > 12.
## Use print(x, correlation=TRUE) or
## vcov(x) if you need itdrop1(CEWL_mod2)## Single term deletions
##
## Model:
## TEWL_g_m2h ~ humidity_tmt_percent * (day + region) + cloacal_temp_C +
## (1 | trial_number/individual_ID)
## npar AIC
## <none> 2563.5
## cloacal_temp_C 1 2619.8
## humidity_tmt_percent:day 1 2622.6
## humidity_tmt_percent:region 4 2558.4We can drop the humidity:region interaction.
CEWL_mod3 <- CEWL %>%
dplyr::filter(complete.cases(.)) %>%
lmerTest::lmer(data = .,
TEWL_g_m2h ~
day*humidity_tmt_percent + region + cloacal_temp_C +
(1|trial_number/individual_ID))
summary(CEWL_mod3)## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: TEWL_g_m2h ~ day * humidity_tmt_percent + region + cloacal_temp_C +
## (1 | trial_number/individual_ID)
## Data: .
##
## REML criterion at convergence: 2510.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3478 -0.6058 -0.1117 0.4446 3.9319
##
## Random effects:
## Groups Name Variance Std.Dev.
## individual_ID:trial_number (Intercept) 38.023 6.166
## trial_number (Intercept) 8.872 2.979
## Residual 123.111 11.096
## Number of obs: 326, groups: individual_ID:trial_number, 33; trial_number, 4
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -69.3758 11.8197 230.9696 -5.869 1.51e-08
## dayAfter 18.2632 1.9251 304.9541 9.487 < 2e-16
## humidity_tmt_percentDry 0.7605 2.7927 41.1331 0.272 0.7867
## regionVentrum 15.7651 1.9396 286.0095 8.128 1.33e-14
## regionHead 2.9138 1.9396 286.0095 1.502 0.1341
## regionDewlap 3.0977 1.9483 286.0970 1.590 0.1130
## regionMite Patch 3.7912 1.9565 286.2700 1.938 0.0536
## cloacal_temp_C 3.8920 0.4886 279.6344 7.966 4.17e-14
## dayAfter:humidity_tmt_percentDry -20.1375 2.4881 288.3429 -8.094 1.63e-14
##
## (Intercept) ***
## dayAfter ***
## humidity_tmt_percentDry
## regionVentrum ***
## regionHead
## regionDewlap
## regionMite Patch .
## cloacal_temp_C ***
## dayAfter:humidity_tmt_percentDry ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) dyAftr hmd__D rgnVnt regnHd rgnDwl rgnMtP clc__C
## dayAfter -0.507
## hmdty_tmt_D 0.005 0.217
## regionVntrm -0.086 -0.004 0.000
## regionHead -0.086 -0.004 0.000 0.504
## regionDewlp -0.098 -0.010 -0.006 0.502 0.502
## reginMtPtch -0.108 0.013 0.002 0.500 0.500 0.498
## clocl_tmp_C -0.972 0.456 -0.122 0.004 0.004 0.017 0.026
## dyAftr:h__D 0.191 -0.680 -0.418 0.004 0.004 0.017 -0.010 -0.146drop1(CEWL_mod3)## Single term deletions using Satterthwaite's method:
##
## Model:
## TEWL_g_m2h ~ day * humidity_tmt_percent + region + cloacal_temp_C + (1 | trial_number/individual_ID)
## Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
## region 9880.0 2470.0 4 286.10 20.063 1.409e-14 ***
## cloacal_temp_C 7812.3 7812.3 1 279.63 63.458 4.165e-14 ***
## day:humidity_tmt_percent 8064.6 8064.6 1 288.34 65.507 1.634e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The model is best with all the parameters currently included in model 3.
write.csv(broom::tidy(CEWL_mod3),
"./best models/exp_effects_CEWL.csv")Change in CEWL
CEWL_before <- CEWL %>%
dplyr::filter(day == "before") %>%
dplyr::select(CEWL_before = TEWL_g_m2h,
humidity_tmt_percent, trial_number,
individual_ID, region)
CEWL_after <- CEWL %>%
dplyr::filter(day == "after") %>%
dplyr::select(CEWL_after = TEWL_g_m2h,
humidity_tmt_percent, trial_number,
individual_ID, region)
CEWL_diffs <- CEWL_before %>%
left_join(CEWL_after, by = c('individual_ID', 'region',
'humidity_tmt_percent', 'trial_number')) %>%
mutate(CEWL_diff = CEWL_after - CEWL_before)plot:
CEWL_diffs %>%
ggplot(data = .) +
geom_boxplot(aes(x = humidity_tmt_percent,
y = CEWL_diff,
group = humidity_tmt_percent,
color = humidity_tmt_percent
),
size = 1,
alpha = 1) +
#facet_wrap(~humidity_tmt_percent) +
theme_classic() +
geom_hline(yintercept = 0, lty = 2) +
xlab("") +
ylab("CEWL Change (g/m2h)") +
#annotate("text", x = 1.5, y = 45,
# label = "paste(italic(p), \" = 0.0152\")",
# parse = TRUE,
# size = 6) +
#ylim(10, 50) +
#scale_x_discrete(labels = c("F" = "Female",
# "M" = "Male")) +
scale_color_brewer(palette = "Set2") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 14),
legend.text.align = 0,
legend.position = "none"
)
Multi-Figure
ggarrange(tmt_effects_osml, tmt_effects_SMI, tmt_effects_hct,
ncol = 1, nrow = 3,
common.legend = TRUE,
legend = "bottom"
) -> tmt_multi_fig## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).## Warning: Removed 3 rows containing missing values (geom_point).## Warning: Removed 3 row(s) containing missing values (geom_path).tmt_multi_fig
# export figure
ggsave(filename = "tmt_multi_fig.jpeg",
plot = tmt_multi_fig,
path = "./final_figures",
device = "jpeg",
dpi = 1200,
width = 6, height = 12)Rehydration…
Data
First, get only the data for before experiment, after experiment, and after rehab.
summary(all_dat)## date individual_ID mass_g hematocrit_percent
## Min. :2021-04-19 37 : 6 Min. : 6.70 Min. :12.00
## 1st Qu.:2021-04-30 39 : 6 1st Qu.:10.20 1st Qu.:24.00
## Median :2021-05-07 40 : 6 Median :11.50 Median :30.00
## Mean :2021-05-06 49 : 6 Mean :11.27 Mean :29.58
## 3rd Qu.:2021-05-13 52 : 6 3rd Qu.:12.60 3rd Qu.:35.00
## Max. :2021-05-20 47 : 5 Max. :15.00 Max. :54.00
## (Other):116 NA's :12
## type osmolality_mmol_kg temp_tmt_C humidity_tmt_percent trial_number
## exp :82 Min. :298.0 25:151 Humid:77 1:35
## rehab :34 1st Qu.:342.8 Dry :74 2:24
## capture:35 Median :359.0 3:44
## Mean :362.5 4:48
## 3rd Qu.:379.0
## Max. :441.0
## NA's :15
## conclusion SVL_mm capture_date day
## complete:151 Min. :59.00 Min. :2021-04-19 Min. : 0.000
## 1st Qu.:66.00 1st Qu.:2021-04-26 1st Qu.: 2.000
## Median :68.00 Median :2021-05-03 Median : 4.000
## Mean :67.45 Mean :2021-04-30 Mean : 5.424
## 3rd Qu.:70.00 3rd Qu.:2021-05-10 3rd Qu.: 9.000
## Max. :73.00 Max. :2021-05-10 Max. :11.000
##
## SMI
## Min. : 7.343
## 1st Qu.: 8.990
## Median :10.011
## Mean : 9.983
## 3rd Qu.:10.751
## Max. :13.970
## rehydrat_dat <- all_dat %>%
dplyr::filter(day %in% c(0, 8, 9, 10, 11))
rehydrat_dat$day <- factor(rehydrat_dat$day,
levels = c(0, 8, 9, 10, 11),
labels = c("Before Experiment",
"After Experiment",
"After Experiment",
"After Rehydration",
"After Rehydration"))
summary(rehydrat_dat)## date individual_ID mass_g hematocrit_percent
## Min. :2021-04-19 37 : 3 Min. : 6.70 Min. :14.00
## 1st Qu.:2021-05-03 39 : 3 1st Qu.:10.20 1st Qu.:23.88
## Median :2021-05-10 40 : 3 Median :11.40 Median :29.50
## Mean :2021-05-07 49 : 3 Mean :11.35 Mean :29.84
## 3rd Qu.:2021-05-13 52 : 3 3rd Qu.:12.80 3rd Qu.:35.00
## Max. :2021-05-20 54 : 3 Max. :15.00 Max. :54.00
## (Other):86 NA's :12
## type osmolality_mmol_kg temp_tmt_C humidity_tmt_percent trial_number
## exp :35 Min. :298.0 25:104 Humid:53 1:17
## rehab :34 1st Qu.:347.0 Dry :51 2:18
## capture:35 Median :369.0 3:33
## Mean :368.5 4:36
## 3rd Qu.:389.0
## Max. :441.0
## NA's :15
## conclusion SVL_mm capture_date day
## complete:104 Min. :59.00 Min. :2021-04-19 Before Experiment:35
## 1st Qu.:66.00 1st Qu.:2021-04-26 After Experiment :35
## Median :68.00 Median :2021-05-03 After Rehydration:34
## Mean :67.68 Mean :2021-05-01
## 3rd Qu.:70.00 3rd Qu.:2021-05-10
## Max. :73.00 Max. :2021-05-10
##
## SMI
## Min. : 7.343
## 1st Qu.: 8.965
## Median :10.011
## Mean : 9.947
## 3rd Qu.:10.731
## Max. :13.970
## prepare summary data to overlay on plots:
hydrat_means <- rehydrat_dat %>%
dplyr::filter(complete.cases(osmolality_mmol_kg, SMI,
hematocrit_percent)) %>%
group_by(humidity_tmt_percent, day) %>%
summarise(osml_mean = mean(osmolality_mmol_kg),
SMI_mean = mean(SMI),
hct_mean = mean(hematocrit_percent))## `summarise()` regrouping output by 'humidity_tmt_percent' (override with `.groups` argument)hydrat_means## # A tibble: 6 x 5
## # Groups: humidity_tmt_percent [2]
## humidity_tmt_percent day osml_mean SMI_mean hct_mean
## <fct> <fct> <dbl> <dbl> <dbl>
## 1 Humid Before Experiment 375. 10.8 36.4
## 2 Humid After Experiment 362. 9.94 27.9
## 3 Humid After Rehydration 367. 9.72 23.5
## 4 Dry Before Experiment 365. 10.5 35.5
## 5 Dry After Experiment 371. 9.21 26.7
## 6 Dry After Rehydration 371. 9.23 22.7SMI
ggplot() +
geom_point(data = rehydrat_dat,
aes(x = day,
y = SMI,
color = humidity_tmt_percent
),
size = 1,
alpha = 0.6) +
geom_line(data = rehydrat_dat,
aes(x = day,
y = SMI,
group = individual_ID,
color = humidity_tmt_percent),
alpha = 0.2) +
geom_line(data = hydrat_means, # new data
aes(x = day,
y = SMI_mean,
group = humidity_tmt_percent,
color = humidity_tmt_percent),
size = 1.6,
alpha = 1) +
theme_classic() +
scale_color_brewer(palette = "Set2",
name = "Humidity Treatment") +
scale_x_discrete(labels = c("Before\nExperiment",
"After\nExperiment",
"After\nRehydration")) +
xlab("") +
xlab("") +
ylab("Body Condition (g)") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text = element_text(color = "black",
family = "sans",
size = 14),
legend.text.align = 0,
legend.position = "bottom"
) -> rehab_SMI_fig
rehab_SMI_fig
# export figure
#ggsave(filename = "rehab_SMI_fig.jpeg",
# plot = rehab_SMI_fig,
# path = "./final_figures",
# device = "jpeg",
# dpi = 1200,
# width = 5, height = 4)Osmolality
first, make a list of all the IDs that have a post-rehab osmolality measurement, since this has a lot of missing data
rehab_osmols <- rehydrat_dat %>%
dplyr::filter(day == "After Rehydration") %>%
dplyr::filter(complete.cases(osmolality_mmol_kg))rehydrat_dat %>%
dplyr::filter(individual_ID %in% rehab_osmols$individual_ID) %>%
ggplot(data = .) +
geom_point(aes(x = day,
y = osmolality_mmol_kg,
color = humidity_tmt_percent
),
size = 1,
alpha = 0.6) +
geom_line(aes(x = day,
y = osmolality_mmol_kg,
group = individual_ID,
color = humidity_tmt_percent),
alpha = 0.2) +
geom_line(data = hydrat_means, # new data
aes(x = day,
y = osml_mean,
group = humidity_tmt_percent,
color = humidity_tmt_percent),
size = 1.6,
alpha = 1) +
theme_classic() +
scale_color_brewer(palette = "Set2",
name = "Humidity Treatment") +
scale_x_discrete(labels = c("Before\nExperiment",
"After\nExperiment",
"After\nRehydration")) +
xlab("") +
xlab("") +
ylab("Osmolality (mmol/kg)") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text = element_text(color = "black",
family = "sans",
size = 14),
legend.text.align = 0,
legend.position = "bottom"
) -> rehab_osml_fig
rehab_osml_fig
# export figure
#ggsave(filename = "rehab_osml_fig.jpeg",
# plot = rehab_osml_fig,
# path = "./final_figures",
# device = "jpeg",
# dpi = 1200,
# width = 5, height = 4)Hematocrit
first, make a list of all the IDs that have all three measurements:
rehab_hct <- rehydrat_dat %>%
dplyr::filter(day == "After Rehydration") %>%
dplyr::filter(complete.cases(hematocrit_percent))rehydrat_dat %>%
dplyr::filter(individual_ID %in% rehab_osmols$individual_ID) %>%
ggplot(data = .) +
geom_point(aes(x = day,
y = hematocrit_percent,
color = humidity_tmt_percent
),
size = 1,
alpha = 0.6) +
geom_line(aes(x = day,
y = hematocrit_percent,
group = individual_ID,
color = humidity_tmt_percent),
alpha = 0.2) +
geom_line(data = hydrat_means, # new data
aes(x = day,
y = hct_mean,
group = humidity_tmt_percent,
color = humidity_tmt_percent),
size = 1.6,
alpha = 1) +
theme_classic() +
scale_color_brewer(palette = "Set2",
name = "Humidity Treatment") +
scale_x_discrete(labels = c("Before\nExperiment",
"After\nExperiment",
"After\nRehydration")) +
xlab("") +
ylab("Hematocrit (%)") +
theme(text = element_text(color = "black",
family = "sans",
size = 18),
axis.text = element_text(color = "black",
family = "sans",
size = 12),
legend.text = element_text(color = "black",
family = "sans",
size = 14),
legend.text.align = 0,
legend.position = "bottom"
) -> rehab_hct_fig
rehab_hct_fig
# export figure
#ggsave(filename = "rehab_hct_fig.jpeg",
# plot = rehab_hct_fig,
# path = "./final_figures",
# device = "jpeg",
# dpi = 1200,
# width = 5, height = 4)Multi-Figure
ggarrange(rehab_osml_fig, rehab_SMI_fig, rehab_hct_fig,
ncol = 1, nrow = 3,
common.legend = TRUE,
legend = "bottom"
) -> rehab_multi_fig
rehab_multi_fig
# export figure
#ggsave(filename = "rehab_multi_fig.jpeg",
# plot = rehab_multi_fig,
# path = "./final_figures",
# device = "jpeg",
# dpi = 1200,
# width = 6, height = 12)